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In past works, various schemes for adaptive synchronization of chaotic systems have been pro- 

(N 

| posed. The stability of such schemes is central to their utilization. As an example addressing this 

^3 ' issue, we consider a recently proposed adaptive scheme for maintaining the synchronized state of 



o 

(N 



identical coupled chaotic systems in the presence of a priori unknown slow temporal drift in the 
couplings. For this illustrative example, we develop an extension of the master stability function 
technique to study synchronization stability with adaptive coupling. Using this formulation, we 
examine local stability of synchronization for typical chaotic orbits and for unstable periodic orbits 
within the synchronized chaotic attractor (bubbling). Numerical experiments illustrating the results 
are presented. We observe that the stable range of synchronism can be sensitively dependent on the 



adaption parameters, and we discuss the strong implication of bubbling for practically achievable 



O 
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adaptive synchronization. 



■ We consider an adaptive scheme for maintaining the synchronized state in a network of identical cou- 

m 

pled chaotic systems in the presence of a priori unknown slow temporal drift in the couplings. Stability 

OO ■ 

of this scheme is addressed through an extension of the master stability function technique to include 

o : 

adaptation. We observe that noise and/or slight nonidenticality between the coupled systems can 
be responsible for the occurrence of intermittent bursts of large desynchronization events (bubbling). 
Moreover, our numerical computations show that, for our adaptive synchronization scheme, the pa- 



rameter space region corresponding to bubbling can be rather substantial. This observation becomes 
important to experimental realizations of adaptive synchronization, in which small mismatches in the 
parameters and noise cannot be avoided. We also find that, for our coupled systems with adaptation, 
bubbling can be caused by a slow drift in the coupling strength. 



I. INTRODUCTION 



QflO 



It has been shown [1|, |2j, |3J] that, in spite of their random-like behavior, the states Xi(t) (i = 1,2, ...,N) of a 
collection of N interacting chaotic systems that are identical can synchronize (i.e., be attracted toward a common 
chaotic evolution, xi(t) — X2(t) = ... = X]y(t)) provided that they are properly coupled. This phenomenon has been 



the basis for proposals for secure communication 4j, |5|, 6 L system identification 



sensors 



13(, information encoding and transmission 



14 



15j . multiplexing 



1 01 ] , data assimilation 



11 



16J, combatting channel distortion 



etc. In all of these applications it is typically assumed that one has accurate knowledge of the interaction between the 

systems, allowing one to choose the appropriate coupling protocol at each node (here we use the network terminology, 

referring to the N chaotic systems as N nodes of a connected network whose links correspond to the input that 

node i receives from node j). In a recent paper jisj l. an adaptive strategy was proposed for maintaining synchronization 

between identical coupled chaotic dynamical systems in the presence of a priori unknown, slowly time varying coupling 

strengths (e.g., as might arise from temporal drift of environmental parameters). This strategy was successfully tested 

on computer simulated networks of many coupled dynamical systems in which, at each time, every node receives only 

one aggregate signal representing the superposition of signals transmitted to it from the other network nodes. In 

addition, the strategy has also been successfully implemented in an experiment on coupled optoelectronic feedback 
□ 

loop s [191 ]. Furthermore, a more generalized adaptive strategy, suitable for sensor applications, has also been proposed 
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In past works, various other schemes for adaptive synchronization of chaos have also been proposed 
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291 ]. So far, in all these studies, when the question of stability of the considered adap tive schemes 



22 



23 



25 



has been studied, the question has been addressed using the Lyapunov function method (see e.g. 
which provides a sufficient but not necessary condition for stability. While this technique has the advantage that 
it can sometime yield global stability conditions, it also has the disadvantages that its applicability is limited to 
special cases, and its implementation, when possible, requires nontrivial system specific analysis. In this paper, we 
address the stability of adaptive synchronization for the example of the scheme discussed in Ref. 18]. In particular, 
our analysis will extend the previously developed stability analysis of chaos synchronization by the master stability 
function technique [ll Q] to include adaptation. We will observe that the range in which the network eigenvalues are 
associated with stability, is dependent on the choice of the parameters of the adaptive strategy. The type of analysis 
we present, while for a specific illustrative adaptive scheme, can be readily applied to other adaptive schemes (e.g., 



those in 



28, 



As compared to the Lyapunov technique, master stability techniques are much more generally applicable but they 
provide conditions for local, rather than global stability. We also note that, within that context, the master stability 
technique allows one to distinguish between stability of typical chaotic orbits and stability of atypical orbits within 
the synchronizing chaotic attractor [i.e., stability to 'bubbling' 



30 



31 



32 



33, 
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; see Sees. Ill and IV]. 



181 ]. which applies to a network of 



In Sec. II we review the adaptive synchronization strategy formulation of Ref. 
chaotic systems with unknown temporal drifts of the couplings. In Sec. Ill, we present a master stability function 
approach to study linear stability of the synchronized solution in the presence of adaptation; we also consider a 
generalized formulation of our adaptive strategy and study its stability. Numerical simulations are finally presented 
in Sec. IV. Our work in Sec. IV highlights the important effect of bubbling in the dynamics. 



II. ADAPTIVE STRATEGY FORMULATION 



As our example of the application of the master stability technique to an adaptive scheme, we consider the particular 



scheme presented in Ref. 18| . To provide background, in this section we present a brief exposition of a formulation 
similar to that in Ref. as motivated by the situation where the couplings are unknown and drift with time. We 
consider a situation where the dynamics at each of the network nodes is described by, 



±i{t) = Fix.it)) + jT[ai{t)ri(t) - H( Xi {t))], 



> 1 V. 



(1) 



where, a;, is the m-dimensional state of system i = 1, N; Fix) determines the dynamics of an uncoupled (7 — > 0) 
system (hereafter assumed chaotic), F : R m — > R m ; Hix) is a scalar output function, H : R m — » R. We take T to be 
a constant m- vector, T = [Ti, T%, T m ] T , with ^\ rf = 1, and the scalar 7 is a constant characterizing the strength 
of the coupling. The scalar signal each node i receives from the other nodes in the network is, 



n(t)=Y,Mi(t)H{xi(t)). 



(2) 



The quantity Aij(t) is an adjacency matrix whose value specifies the strength of the coupling from node j to node i. 
We note that if 



(3) 
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then Eq. JT]) admits a synchronized solution, 

xi(t) = x 2 (t) = ... = x N (t) = x s (t), (4) 

where x s (t) satisfies 

x s (t) = F(x s (t)), (5) 



which corresponds to the dynamics of an isolated system. We regard the Aij(t) as unknown at each node i, while 
the only external information available at node i is its received signal The goal of the adaptive strategy is to 
adjust (Ji(t) so as to maintain synchronism in the presence of slow, a priori unknown time variations of the quantities 
A.y(t). That is, we wish to maintain approximate satisfaction of Eq. ([3]). For this purpose, as discussed in Ref. 18J, 
our scheme can be extended to the case where the output function is ^-dimensional, H : R m — > R e , where £ < m 
and r is an £ x m dimensional matrix. For simplicity we consider £ = 1. We assume that each node independently 
implements an adaptive strategy. At each system node z, we define the exponentially weighted synchronization error 
ipi =< (a t ri — H{xi)) 2 >„, where 



< 



G(t) > u = / Gity-^'-^dt', (6) 



and we evolve <7i(t) so as to minimize this error (a slightly more general approach is taken in 18|). Hence we set 
dipi/dai equal to zero to obtain, 

< H(x t (t)) ri (t) > v Pi (t) 

By virtue of d < G(t) > l/ /dt = — v < G(t) -\-G(t), we obtain the numerator and the denominator on the right 
hand side of Eq. ([7]) by solving the differential equations, 

Pi(t) = -vpi(t) + n(t)H(xi(t)), (8a) 
ft(t) = -v qi {t) + ri {t) 2 . (8b) 



Since the dynamics of A^ (t) is imagined to occur on a timescale which is slow compared to the other dynamics in the 
network, we can approximate (t) as constant Ay . This essentially assumes that we are dealing with perturbations 
from synchronization whose growth rates (in the case of unstable synchronization) or damping rates (in the case of 
stable synchronization) have magnitudes that substantially exceed \A~^(t)dAi 3 ■ /dt\. Under this assumption, we note 
that Eqs. {!]), ©i an d ® admit a synchronized solution, given by Eqs. J4}, ©, and 



Pi = -"Pi + (J2^)H(x s ) 2 , i = 1, ...,N, 

3 

Qt = '•'/; + (£A ij )*H{x')*, i = 1 V. 



(9a) 
(9b) 



To simplify the notation, in what follows, we take DF s (t) = DF(x s (tj), H s {t) = H{x s (t)), and DH s (t) = DH(x s (t)); 
e.g., we can now write, 



v\ = k< {H s f > v , 
q? = k\ < (H s f > v , 



(10) 



where fe« = {^2,- Aij). If the synchronization scheme is locally stable, we expect that the synchronized solution (|3]),([5]), 
and ^ will be maintained under slow time evolution of the couplings A^ (t) . 

III. STABILITY ANALYSIS 
A. Linearization and master stability function 

Our goal is to study the stability of the reference solution ([I]),©, and By linearizing Eqs. |T]) and © about 
(O, and we obtain, 



S±i = DF s Sx t +-/Tl DH 



h = -vet - H s DH s kj 



i=l, ...,7V, 
i=l,...,N, 



(11a) 

(lib) 



where we have introduced the new variable ej(i) = kiSpi(t) — Sqi(t). 



Equations (fTTjl constitute a system of (m + l)N coupled equations. In order to simplify the analysis, we seek to 
decouple this system into N independent systems, each of dimension (m + 1). For this purpose we seek a solution 
where 5xt is in the form 5xi = Cix(t), where Ci is a time independent scalar that depends on i and x{t) is a m-vector 
that depends on time but not on i. Substituting in Eqs. (|lla[) . (|llb[) . we obtain, 



x = DF s x + 7 r 



- 1 



DH s x 



~fTH s 



i = l,...,N, 



-va - ki 



^ ^ Aij Cj ki Ci 



H s DH s x, i = l,...,N. 



(12a) 
(12b) 



To make Eqs. (1121) independent of i, we consider (3(t) — ei(t)/[ciki (a — 1)] and Yl- AtfCj = akiCt, where a is a 
quantity independent of i. Namely, the possible values of a are the eigenvalues, A'c = qc, corresponding to linearly 
independent eigenvectors c = [ci,C2, ...,cjv] t , where A' = {A'^} — {k i ~ 1 A i j}. This gives, 



x = DF s x - 7(1 - a) 



TDH s x + T 



H s f3 



v(3 - H s DH s x 



(13a) 
(13b) 



which is independent of i, but depends on the eigenvalue a. Considering the typical case where there are N distinct 
eigenvalues of the NxN matrix A', we see that Eqs. (fT2)) constitute N decoupled linear ordinary differential equations 
for the synchronization perturbation variables x and f3. All the rows of A' sum to 1. Therefore A' has at least one 
eigenvalue a — 1, corresponding to the eigenvector c\ = C2 = ... = cjy = 1- Furthermore, since A'^ > for all 
we have by the Perron-Frobenius theorem that a < 1, and thus (1 — a) > 0. For a = 1, Eq. (|13ap becomes, 



x = DF s x. 



(14) 



This equation reflects the chaos of the reference synchronized state (Eq. ([5])) and (because all the a are equal) 
is associated with perturbations which are tangent to the synchronization manifold and are therefore irrelevant 
in determining synchronization stability. Stability of the synchronized state thus demands that Eqs. (|12|) yield 
exponential decay of x and (3 for all the (N — 1) eigenvalues a, excluding this a = 1 eigenvalue. 

Then it becomes possible to introduce a master stability function l|, [3j, M(£), that associates the maximum 
Lyapunov exponent of system (fT3"|) with £ = 7(1 — a). In so doing, one decouples the effects of the network topology 



(reflected in the eigenvalues a and hence the relevant values of £ = 7(1 — a)) from the choices of F, H, v. In general 
an eigenvalue, and hence also £, can be complex. For simplicity, in our discussion and numerical examples to follow, 
we assume that the eigenvalues are real (which is for instance the case when the adjacency matrix is symmetric). 
For any given value of 7 stability demands that M(£) < for all those values of £ = 7(1 — a) corresponding to the 
eigenvalues a 1. 
Following Refs. 



36 



37 



38 



39|, we now introduce the following definition of synchronizability. Let us assume that 
the master stability function M(£) is negative in a bounded interval of values of £, say [£~, £ + ]. Then, in order for the 
network to synchronize, two conditions need to be satisfied, (i) £~ < 7(1 — a m i n ), and (ii) £ + > 7(1 — a ma x), where 
dmin {czmax) is the smallest (largest) network eigenvalue over all the eigenvalues 1. The network synchronizability 
is defined as the width of the range of values of 7, for which Af (£) < 0. Assuming that a m i n and a max are assigned 
(e.g., the network topology is given), then the network synchronizability increases with the ratio £ + /£~. In what 
follows, we will compare different adaptive strategies in terms of their effects on the synchronizability ratio £ + /£~- 

In our analysis above, since we divide by fej, we have implicitly assumed that all the fej 7^ 0, i.e., that every node 
has an input. There is, however, a case of interest where this is not so, and this case requires separate consideration. 
In particular, say there is one and only one special node (which we refer to as the maestro or sender) that has no 
inputs, but sends its output to other nodes (which interact with each other), and we give this special node the label 
i = N. Since node N receives no inputs, we do not include adaption on this node, and we replace Eq.(l) for i = N by 
iivCO = F(%N(t)). In addition, when investigating the stability of the synchronized state, it suffices to set Sx^it) = 
(i.e., not to perturb the maestro). Following the steps of our previous stability analysis, we again obtain Eqs. (|12p 
and (|13[) . but with important differences. Namely, Eqs. (|12p now apply for i = 1, ..., N — 1, the values of a in Eqs. 
(fT3|) are now the eigenvalues of the (N — 1) x (JV — 1) matrix {A^} = {jfej Ay} for i, j — 1, 2, (JV — 1); i.e., only 
the interactions between the nodes i,j < (N — 1) are included in this matrix. Note that fcj is still given by X)j=i -^»i> 
still including the input Ajjv from the maestro node. Also since 5xn = 0, all of the eigenvalues represent transverse 
perturbations and are therefore relevant to stability. (This is in contrast to the case without a maestro in which we 
had to exclude an eigenvalue, i.e., a = 1 corresponding to ci = C2 = ... = cn = 1. For a similar discussion for the 



40(.) The simplest case of this type (used in 



case of the standard master stability problem with no adaptation, see 
some of our subsequent numerical experiments) is the case N = 2, where there is one receiver node (i = 1) and one 
sender/maestro node (i = 2). Since there is only one receiver node whose only input is received by the sender, A 
reduces to the scalar A — and a = 0, yielding £ = 7. 



8 

As stated above, x s (t) in ((4]) is an orbit of the uncoupled system (|4]). In general, two types of orbits x s (t) are of 
interest: (i) a typical chaotic orbit on the relevant chaotic attractor of (j4|), and (ii) the orbit that is ergodic on the 
maximally synchronization- unstable invariant subset embedded within the relevant chaotic attractor of (|4]). Here, by 
'relevant chaotic attractor' we mean that, if the system (QJ has more than one attractor, then we restrict attention to 
that attractor on which synchronized motion is of interest. Also, in (i), by the word 'typical', we mean orbits of (j4]) 
that ergodicly generate the measure that applies for Lebesgue almost every initial condition in the attractor's basin 
of attraction. In this sense, the orbit in (ii) is not typical. In general the criterion for stability as assessed by (ii) 
is more restrictive than that assessed by (i). Conditions in which the synchronized dynamics is stable according to 



(i), but unstable according to (ii), are referred to as the 'bubbling' regime 



synchronization of chaos 



30 



31 
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35|| . In previous work on 



35j, it has been shown that, when the system is in the bubbling regime, 
small noise and/or small 'mismatch' between the coupled systems can lead to rare, intermittent, large deviations from 
synchronism, called 'desynchronization bursts' 



42j. By small system mismatch we mean that, for each node i, the 
functions F in {1} are actually different, F — » Fi, but that these differences are small (i.e., \Fi(x) — F(x)\ is small, 
where F(x) now denotes a reference uncoupled system dynamics; e.g., Fi averaged over i). With reference to our 
adaptive synchronization problem |T]), we shall see that, in addition to small noise and small mismatch in F, bursting 
can also be induced by slow drift in the unknown couplings Aij(t). From the practical, numerical perspective, the 
complete and rigorous application of the stability criterion (ii) is impossible, since there will typically be an infinite 
number of distinct invariant sets embedded in a chaotic attractor, and, to truly be sure of stability, each of these 
must be found and numerically tested. In practice, therefore, as done previously by others, we will evaluate stability 
for all the unstable periodic orbits embedded in the attractor up to some specified period. This will give a necessary 



33] 



condition for stability according to (ii), and furthermore, it has been argued and numerically verified in Ref. 
that stability, as assessed from a large collection of low period periodic orbits (and embedded unstable fixed points, if 
they exist in the relevant attractor), will extremely often yield the true delineation of the parameters of the bubbling 
regime, or, if not, an accurate approximation of it. Our numerical results of Sec. IV lend further support to this idea. 



B. Generalized adaptive strategy 



We now analyze a generalization of our adaptive strategy. Namely, we replace Eq. (|8b[) by. 



nit) = -m(t) + h(t)/p t (t)]H(x t (t)YQ 



2n / Pi(t)ri(t) 



q % {t)H{ Xl {t)) 



(15) 



where Q(z) is an arbitrary function of z, normalized so that Q(l) = 1. The key point is that at synchronism 
OiTi = H(xi(t)), corresponding to p^n = qiH(xi(t)); and thus, since we take Q(l) = 1, the synchronized solution is 
unchanged. The stability analysis for this generalization is given in the Appendix I and results in the following master 
stability equations, 



x = DF s x ~ £ 



ITS Of 

TDH s x + T- 



< (# s ) 2 >, 



(16a) 



/3> = -v0 + (0 - 1) ( ^ )2 + (0 - 2)H s DH s x, (16b) 

where = Q'(l), and Q'(l) denotes dQ(z)/dz evaluated at z = 1. We then introduce a master stability function 
M(£,<j>), that associates the maximum Lyapunov exponent of system (|16|) with £ = 7(1 — a) and <f>. 

Thus we expect that, when our modified adaptive scheme is stable, it will again relax to the desired synchronous 
solution. The difference between the stability of the modified scheme (Eqs. (TT5)) ') and the stability of the original 
scheme (corresponding to Eqs. (|16p with <p — 1), is that, by allowing the freedom to choose the value of <j), we can 
alter the stability properties of the synchronous state. We anticipate that, by properly adjusting </>, we may be able 
to tailor the stability range to better suit a given situation. 

In the case of (f> = 2, Eq. (|16b[) reduces to, 



which has a Lyapunov exponent A = Aq — v, where Ao is the time average of (H s ) 2 / < (H s ) 2 >„, Ao > 0. For v > Ao, 
Eq. (jXTJ) implies that /?' decays to zero. Thus, if we choose a large enough value of 1/, stability of the synchronized 
state is determined by (|16ap with (3 1 set equal to zero, and Eq. (116a[) reduces to the master stability function for the 
determination of the stability of the system without adaptation [3|. Therefore, in the case of <fi = 2, v > Ao, the stable 
range of 7 is independent of v and is the same as that obtained for the case in which adaptation is not implemented 
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IV. NUMERICAL EXPERIMENTS 

In our numerical experiments we consider the example of the following Rossler equation, for which, m — 3, x(t) 
(u(t),v{t),w(t)) T , 



F(x) = 



-v — w 



u + av 
b + (u — c)w 



(18) 



with the parameters a = b = 0.2, and c = 7, and we use H(x(t)) = u(t), and T = [1,0, 0] T . In Fig. 1 the master 
stability functions M(£) calculated from Eq. (fT3"|) for the adaption scheme of Sec. II are plotted for three different 
values of is, i.e., v = 0.1,2,6 (dashed, dashed/dotted, and dotted curves, respectively). In addition, for comparison, 
we also plot the result of M(£) computations for the case in which no adaptation is introduced, corresponding to 
the reduced system x = [DF S + j(a — l)TDH s ]x (solid curves). The master stability function is shown in black 
(respectively, grey) for the cases that x s (t) is a typical chaotic orbit in the attractor (respectively, the maximally 
unstable periodic orbit embedded in the attractor for periodic orbits of period up to four surface of section piercings; 
see Appendix II for a brief account of how the unstable periodic orbits were obtained) . We say that synchronization 
is 'high quality' stable in the range of £ for which M(£) for all orbits (i.e., including the periodic orbits) is negative. 
As can be seen, by changing the parameter is, the £-range of stability can be dramatically modified. The bubbling 
range is given by the values of £ for which M(£) < for a typical orbit but M(£) > for the maximally unstable 
periodic orbit embedded in the attractor. 

Figure 2 is a £ — v level curve plot of the values assumed by the master stability function M evaluated for x s (t) 
being a typical chaotic orbit. In the figure, the area of stability (corresponding to M < 0) is delimited by the thick 
0- level contour line. From the figure, we see that the width of the range of stability increases with v. In Figs. 3(a,b) 
a comparison between the areas of stability is given, for the cases in which x s (t) is a typical chaotic orbit in the 
attractor, and for the case that x s (t) is the maximally unstable periodic orbit embedded in the attractor of period up 
to four. The thick solid (respectively, dashed) curves bound the area in which the master stability function M(£, v) is 
negative for x s (t) corresponding to a typical chaotic orbit in the Rosller attractor (respectively, for x s (t) corresponding 
to the maximally unstable periodic orbit embedded in the attractor of period up to four). The bubbling area falls 
between the dashed and the continuous contour lines. 
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FIG. 1: The plot shows the master stability function M(£) versus £ for for the case in which no adaptation was introduced, 
corresponding to a = 1 (black continuous line) and for three different values of u, i.e., v — 0.1, 2, 6 (dashed and dotted lines). 
The master stability functions obtained by choosing x s (t) to be a typical chaotic orbit in the attractor (respectively, the 
maximally unstable periodic orbit embedded in the attractor of period up to four) are in black (respectively, grey). F(x) is the 
Rossler equation ([18]), H(x(t)) = u(t), and T = [1,0, 0] T . 



Interestingly, we see that for 1.2 < v < 3.2, high-quality stability can never be achieved for any £, while, in contrast, 
stability with respect to typical chaotic orbits (i.e., with bubbling) is achievable. Let £ t + , £ t ~~ , , £~ denote the upper 
(+) and lower (— ) values of £ at the borders of the stability regions with respect to a typical (t) chaotic orbit and 
with respect to unstable periodic orbits (p) in the synchronizing attractor. E.g., high-quality synchronism applies for 
> £ > £~ and the bubbling regime corresponds to £~ > £ > £jT or £j > £ > In terms of these quantities, 
useful measures for asses sing the possibility of achieving stable synchronism for a given network topology are the 



'synchronizability' ratios 



3G 



37 



38 



3a, 



£+ 



(19) 



In what follows, where convenient, we drop the subscripts t and p with the understanding that the discussion may 
be taken to apply to stability based on either typical or periodic orbits. Noting that synchronism is stable for 
£ + > £ > £~ , and that £ = 7(1 — a), we consider the coupling network topology-dependent ratio (1 — cv~)/(l — a + ) 
where a + (<x~) denotes the maximum (minimum) eigenvalue of the adjacency matrix (not including the eigenvalue 
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a = 1 corresponding to the eigenvector (1, 1, 1) T ). Recall that (1 — a) > 0. Since £ + > £ > £ for stability, if 

S > jj— £p (20) 

then the system can be made stable by adjustment of the constant 7, but, if s < (1 — a~)/(l — a + ), then it is 
impossible to choose a value of 7 for which M(£) < for all the relevant eigenvalues a, and stability is unachievable. 
Figure 3(c) shows plots of s t and s p versus v for the same parameters as used in Figs. 3(a,b). Note that, for these 
computations, the values of s without adaption (i.e., St = 23.5 and s p = 10.5) always exceed the corresponding values 
with adaption. We have also found this to be true for the generalized adaptive scheme of Sec. Ill B (which includes 
the additional adaption parameter <p). However, we do not know whether this is general, or is limited to our particular 
example (Eq. (|18p with H(x) = u, and our choices of the parameters a,b, and c). 

To test our linear results in Fig. 3, we have also performed fully nonlinear numerical simulations for a simple 
network consisting of a sender system (labeled 1) connected to a receiver (labeled 2). In this case Eq. (1) becomes 

x x {t) =F{x x {t)), (21a) 
i 2 (t) =F(x 2 (t)) + jT[a(t)A(t)H(x 2 (t)) - H( Xl (t))}, (21b) 

and A(t) is a scalar. Each data point shown in Figs. 3(a,b) corresponds to a run, where the sender was given a 
random initial condition and random values for v and £ were chosen in the plotted range. After waiting sufficient 
time to ensure that the sender state is essentially on the attractor, the it-variable of the receiver state was initialized 
by a displacement of 10 -8 from the M-variable of the sender state. A step-size of 10~ 4 was used for a run time of 10 5 
time units, over which we recorded the normalized synchronization error, 

E(t) = (22) 

< {u s — < u s >y > L ' Z 

where < ... > indicates a time average and the subscript s denotes evolution on the synchronous state (i.e., using 
dynamics from Eq. id}). If? hi that time span, E never converged to and, at some point, exceeded 0.1, the run 
was considered to be unstable (corresponding to an x in the hgure). If E converged to 0, a 1% mismatch in the 
Rossler parameter a was introduced to the receiver, and the run of duration 10 5 time units was repeated with an 
initial separation of 0. If, at any time during the run, E ever exceeded 0.1, the run was considered to be bubbling 
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FIG. 2: The figure is a level curve plot in £-1/ space of the values assumed by the master stability function M, evaluated for 
x s (t) being a typical chaotic orbit. The area of stability (corresponding to M < 0) is delimited by the thick 0- level contour 
line. F(x) is the Rossler equation [fig]). H(x(t)) = u(t), and V = [1, 0, 0] T . 

(corresponding to a green circle in the figure), otherwise the run was considered to be stable (corresponding to a 
red triangle in the figure). We see that the master stability computations of the high-quality stable, bubbling, and 
unstable regions (the solid and dashed lines) correspond well with these results. We also did a sampling of points up 
to period 5 and did not find that this altered our results. From Fig. 3(a), we observe the presence of a few green 
circles (i.e, bubbling) within the high-quality synchronization area, delimited by the dashed line. In reference to this 
observation, we note that (i) for the case in which a small parameter mismatch is present, the synchronization error is 
expected to vary smoothly with parameter variation, and there is no sharp transition from the stable to the bubbling 
regime; and (ii) our computations show that close to the dashed line, the master stability function associated with 
the most unstable invariant set embedded in the attractor is rather small. Facts (i) and (ii) explain our difficulty 
in using our nonlinear computations to clearly separate the bubbling from the stable regions about the dashed line 
in Fig. 3(a). An important point concerning Figs. 3(a,b) is that the area associated with bubbling in Fig. 3(a) is 
rather substantial. This observation would become particularly important in experimental realizations of adaptive 
synchronization, since small mismatches in the parameters and noise cannot be avoided in experiments. 

Figure 4 shows a sample plot of the normalized synchronization error E(t), versus t. We implemented our adaptive 
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strategy with values of v — 2.5 and 7 = 5, corresponding to the bubbling regime (see Fig. 3), A(t) = 1, and the 
receiver has 0.1% mismatch in the parameter a. The two insets are zooms showing phase-space projections in the plane 
(112,1)2), over two different time intervals. Inset (b) corresponds to a range of time between bursts (E(t) < 5 x 10 -2 ) 
and shows that during this time the orbit is essentially that of a typical chaotic orbit. Inset (a) shows the orbit 
trajectory for a range of time during which a burst is growing. It is seen from inset (b) that during the time range of 
the growing burst the orbit closely follows a period 4 orbit embedded in the attractor. The burst is evidently caused 
by the instability of this period 4 orbit to perturbations that are transverse to the synchronization manifold. 

We have also performed numerical master stability computations for our generalized adaptive strategy, presented 
in Sec. Illb. This is shown in Fig. 5, where the £ + (</>) an d £ _ (</ ) ) curves, corresponding respectively to the largest 
(smallest) values of £ for which M (£,(/>) > (Af(£, <j>) < 0), are plotted versus </> for three different values of v = 
[0.1,2.0,6.0] for typical chaotic orbits. For small v (e.g., v = 0.1 in the figure), the range of stability is 
almost independent of <fi, while for larger values of v the choice of <fi can significantly affect the £-range of stability. 
As expected, at <j> = 2, £ + (0) and £~(<j>) are independent of v. 

Finally, we investigated whether, for our coupled systems with adaptation, bubbling can be caused by a slow drift 
in the coupling strength. For this purpose we now take the parameter A(t) in Eq. (|2"Tj) to have a slow time drift, 

A(t) = 1 + 0.2 sin(27r x 10~ 3 i). (23) 

We implemented our adaptive strategy with values of v = 1 and 7 = 2, corresponding to the bubbling regime (see 
Fig. 3). For most of the time there is good synchronization between the sender and the receiver, but we also observed 
the intermittent occurrence of short, intense desynchronization bursts. Figure 6 shows the synchronization error 
E[t) versus t. Note that in the absence of parameter drift (A constant), the synchronization error would eventually 
become zero. This simulation shows that, similarly to the previously reported burst-inducing effect of small parameter 
mismatch or noise, drift also promotes the continuous intermittent occurrence of bursting. 

V. CONCLUSION 



This paper is concerned with the study of stability of adaptive synchronization of chaos in coupled complex networks 
(e.g., sensor networks). As an example addressing this issue, we consider a recently proposed adaptive scheme for 
maintaining synchronization in the presence of a priori unknown slow temporal drift in the couplings 18 1 . In contrast 
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with previous approaches (e.g., 



22|, |23j, [25J, |27j ) , based on system specific use of the Lyapunov function technique, we 
present a master stability analysis which predicts the exact ranges of stability for the synchronized state. We observe 
that the stable range of synchronism can be sensitively dependent on the adaption parameters. Moreover, we are able 
to predict the onset of bubbling, which occurs when the synchronized state is stable for typical chaotic orbits but is 
unstable for certain unstable periodic orbits within the synchronized chaotic attractor. We define stability to be high 
quality when the synchronized state is stable with respect to all the orbits embedded in the attractor and numerically 
find the regions of 'high quality stability' for a given system of interest. We also found that, for our coupled systems 
with adaptation, bubbling can be caused by a slow drift in the coupling strength, in addition to small noise and 
small mismatch in F. We emphasize that, since parameter mismatch, noise and drift are ubiquitous in experimental 
situations, and since (e.g., Fig. 3(a)) bubbling can occupy substantial regions of parameter space, consideration of 
bubbling can be expected to be essential for determining the practical feasibility of chaos synchronization applications. 
We thank Anurag. V. Setty and Bhargava Ravoori for the enlightening discussions. 
This work was supported by ONR grant N00014-07- 1-0734. 



APPENDIX I: Stability of the generalized adaptive strategy 

We note that the function Q([pi(t)ri(t)]/[qi(t)H(xi(t))]) in Eq. (IT5|> . when evaluated about (J9j) , is equal to one. 
Then, by linearizing Eqs. ([1]), (|5a|) . and ([T5|) about (O, we obtain, 



Sx t = DF s 8x i +'yT^DH l 



-VCi + (0 - 1) 



(H 



s\2 



< (£P) 2 > v 



2)H S DH S 



kf < (IP) 2 > v 1 



i = l,...,N, 



i = l,...,N, 



(24a) 
(24b) 



As in our derivation of Eqs. (|13p . we again set Sxi = Cix(t), where Ci is a constant scalar that depends on i and x(t) 
is a vector that depends on time but not on i. Equations (|24p . then become 



x = DF s x + 7 r 



DH s x 



jTH s 



-va + {cj> - I) 



s\2 



< (£P) 2 > v 



2)H S DH S 



h% ^ Aij Cj Ci 



ei, i = l,...,N, 



x, 1 = 1, N. 



(25a) 
(25b) 



To make Eqs. ([231) independent of i, we again consider [3'{t) = €i{t) / [ciki(a — 1)] and take a to be the eigenvalues of 
A ' = i A 'ij} = {h^Aij}, resulting in Eqs. ([16]). 
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APPENDIX II: Determination of unstable periodic orbits 

To account for the phenomenon of bubbling, it is necessary to look not just at typical (that is, chaotic) orbits of 
the uncoupled oscillator, but the periodic orbits embedded in the chaotic attractor as well. As there are a (countably) 
infinite number of such orbits, it is impossible to account for them all. However, as shown by Hunt and Ott, the 

n 

optimal periodic orbits of maximal transverse instability tend to be those of low period [33j. Thus, for our analysis, 
it was found to be sufficient to consider only those orbits with a period less than some appropriately chosen limit. 

To find these low-period orbits for the Rossler attractor, we initialized an uncoupled oscillator with random initial 
conditions, waited for it to settle onto the attractor, then recorded its orbits for some suitable length of time at high 
temporal precision. We then noted each piercing of the surface of section u — in the positive- m direction (ii > 0). 
To a high degree of approximation, the (w, w) coordinates of these points were found to lie a curve, thus suggesting 
that it is possible to reduce the three-dimensional flow to a one-dimensional map. We then plotted v(i + n) vs. v(i); 
that is, the v coordinate of the (i + n) th piercing versus the v coordinate of the i th . Each intersection of this curve 
with the line v(i + n) — v(i) represents the v coordinate of an initial condition for an orbit that starts on the surface 
of section and returns to its original position after n piercing of the surface of section. With two coordinates (namely 
u and v) known, all that remains is to find the value of w such that (0, v, w) lies on the attractor. 

Of course, for n > 1, many of these intersections will be redundant, as every period n orbit pierces the surface 
of section n times, thus producing n intersections on the curve. In addition, each curve will have intersections 
corresponding to orbits of any period that is a factor of n. As an example, consider the curve v(i + 4) vs. v(i). The 
Rossler system used in this paper has three Period 4 orbits, one Period 2 orbit and one Period 1 orbit. Thus, the 
number of times v(i + 4) vs. v(i) will intersect v(i + 4) = v(i) is 3x4+1x2 + 1x1 = 15. 

As these orbits are inherently unstable, error accumulated through numerical integration can result in a trajectory 
leaving the periodic orbit after only a small number of periods. Thus, for the long term computation of Lyapunov 
exponents to obtain the master stability function, it is advisable to compute the trajectory for only a single period, 
then return the oscillator to its initial position (it, v,w), and repeat as often as needed. 
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FIG. 3: In plot (a), thick solid curves (thick dashed curves) bound the area in which the master stability function M(£, v) is 
negative for x s (t) corresponding to a typical chaotic orbit in the Rossler attractor (for x s (t) corresponding to the maximally 
unstable periodic orbit embedded in the attractor of period up to four), F(x) is the Rossler equation (|18[) . H(x(t)) = u(t), and 
r = [1, 0, 0] T . Each data point shown in the figure is the result of a simulation involving a sender (maestro) system connected 
to a receiver, where the receiver state was initialized by a displacement of 10 -8 from the sender state. A step-size of 10 -4 was 
used for a run time of 10 time units. If, in that time span, the synchronization error E never converged to and, at some 
point, exceeded 0.1, the run was considered to be unstable (corresponding to an x in the figure). If E converged to 0, a 1% 
mismatch in the Rossler parameter a was introduced to the receiver, and the run was repeated with an initial separation of 0. 
Then, if E ever exceeded 0.1, the run was considered to be bubbling (corresponding to a green circle in the figure), otherwise 
the run was considered to be stable (corresponding to a red triangle in the figure). Plot (b) is a blow up of the lower left corner 
of plot (a). Plot (c) shows the synchronizability ratios st (solid curve) and s p (dashed curve) versus v. The missing data points 
for the dashed curve are a result of the low period orbits not having a range of stability for those values of v. We found that 
the synchronizability ratios for the nonadaptive case are equal to those in the limit v — *• 0. 
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FIG. 4: The figure shows the synchronization error E(t) versus t for a simple network consisting of a sender connected to a 
receiver (Eqs. J2TJ), F(x) is the Rossler equation JJS), H(x(t)) = u(t), T = [1,0, 0] T , 7 = 5, v = 2.5, A(t) = 1, dt = 10 -3 . 
The receiver has a 0.1% mismatch in the parameter a. The two insets are zooms showing phase-space projections in the plane 
(m2,i>2), over two different time intervals. Inset (b) corresponds to a typical chaotic orbit for which the synchronization error 
is small, i.e., E(t) < 5 x 10~ 2 , while inset (a) corresponds to an unstable period 4 periodic orbit embedded in the attractor, for 
which E(t) is eventually large (i.e., a burst occurs). 
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FIG. 5: The plot shows the area in the parameter space (</>, £) in which M(£, </>) obtained from (|16[) is negative, for three 
different values of v = [0.1,2.0,6.0]; Fix) is the Rossler equation QH), H(x(t)) = u(t), and F = [1,0, 0] T . The stability areas 
are upper and lower bounded by the £ + curve and the £~ curve, plotted as function of 0. As the figure shows, at <f> = 2, £ + 
and £~ are independent of z/, corresponding to the case of no-adaptation. 
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FIG. 6: The figure is a plot of the synchronization error E(t) (defined in Eq. (|22[) ) versus t for a simple network consisting of 
a sender connected to a receiver (Eqs. (|21[) ~1. F(x) is the Rossler equation (|18p . H{x{t)) = u(t), F = [1,0, 0] T , u = 1, 7 = 2, 
A(t) = l + 0.2sin(27r x 10~ 3 t), dt = 10 -3 . As can be seen, the dynamics of E(t) exhibits intermittent bursting. 



